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We investigate the effective properties (conductivity, diffusivity and elastic moduli) of model 
random composite media derived from Gaussian random fields and overlapping hollow spheres. The 
morphologies generated in the models exhibit low percolation thresholds and give a realistic rep- 
resentation of the complex microstructure observed in many classes of composites. The statistical 
correlation functions of the models are derived and used to evaluate rigorous bounds on each prop- 
erty. Simulation of the effective conductivity is used to demonstrate the applicability of the bounds. 
The key morphological features which effect composite properties are discussed. 
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I. INTRODUCTION 

The prediction of effective properties of heterogeneous 
systems such as porous media and two phase compos- 
ites is of considerable interest pl-Pl . Understanding the 
inter-relationships between rock properties and their ex- 
pression in geophysical and petrophysical data is neces- 
sary for enhanced characterisation of underground reser- 
voirs. This understanding is crucial to the economics of 
oil and gas recovery, geothermal energy extraction and 
groundwater pollution abatement. Manufactured com- 
posites such as foamed solids H] and polymer blends 0] 
often exhibit a complex microstructure. To optimize the 
properties of these systems it is necessary to understand 
how morphology influences effective properties. In gen- 
eral, the difficulty of accounting for microstructure has 
made exact prediction impossible in all but the simplest 
of cases. 

On the other hand, considerable progress has been 
made in the derivation of rigorous bounds on a host of 
properties HH]. For example, relatively accurate bounds 
have been derived for the elastic moduli and conductiv- 
ity of isotropic two-phase composites ]7|-pd| . To evaluate 
these bounds for a given system it is necessary to know 




the 3-point statistical correlation function 113 . Due to 
the difficulty of measuring this information |l^-[l5||, a 
number of model media have been proposed for which the 
functions can be explicitly evaluated. These include: cel- 
lular jQ, particulate H and periodic ^\ materials (eg. 
Figs. |^(a)&(b)). The principal problem with these mod- 
els is that they employ over-simplified representations of 
the inclusion (or pore) structure observed in many natu- 
ral and manufactured composite materials. 

Recently we derived the properties of a model of amor- 
phous materials |Q (e.g. Fig. |l|(c)) based on level-cut 
Gaussian random fields ||l^ (GRF). Although the GRF 
model is applicable to many classes of non-particulate 
composite materials, it cannot account for materials 
which remain percolative at very low volume fractions. 

Porous rocks |^,Q , polymer blends [g| , solid foams |j] 
and membranes provide examples of systems where a sin- 
gle phase remains connected down to low volume frac- 
tions. The percolation threshold of a system is only 
one factor which determines its effective properties: The 
shape of the pores/inclusions should also be consid- 
ered [pl|-p3[. Polystyrene foam, an example of a highly- 
porous material, is shown in Fig. 0. 





(a) ■■■^'-■- (b) 

FIG. 1. Models of composite microstructure. (a) Periodic models - regular array of spheres; (b) Particulate models - identical 
overlapping spheres (lOS) and; (c) Gaussian random held (GRF) models - single cut variant. None of the models can mimic 
the microstructure observed in percolative low volume fraction materials (e.g. polystyrene foam - Fig. [3). 
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FIG. 2. Morphology of polystyrene foam. 

The complex solid phase has a 'sheet-like' character quite 
different to that found in cellular, particulate |l|| and sin- 
gle level cut GRF |l|,|l|l models (Fig. |i|). It is clear 
that current models of composite microstructure cannot 
account for the percolative and morphological character- 
istics observed in porous rocks, solid foams, membranes 
and polymer blends. 

In this paper we describe models which do give a realis- 
tic representation of the microstructure observed in many 
classes of composite materials, and which remain percola- 
tive at very low volume fractions. Variational bounds and 
computer simulation are used to estimate the influence 
of morphology on diffusive-transport and elastic proper- 
ties. The first model is an extension of the Gaussian 
random field (GRF) model considered in a previous pa- 
per [|l8[ . In this case the interface between the composite 
phases is defined by two (rather than one) level cut of a 
GRF [|9|J2|-|2^. The freedom in choosing the position 
of the cuts (for a given volume fraction), and the spectra 
of the model, allows a rich variety of morphologies to be 
modelled. By qualitatively comparing these morpholo- 
gies to those observed in physical systems the models 
can be associated with classes of physical composites. 

A second highly porous model can be obtained by gen- 
eralizing the well-known lOS model [1| to include the case 
of arbitrarily thin hollow spheres. This model is appli- 
cable to a class of ceramics and foams fabricated from 
hollow spheres: a composite which possesses excellent 
uniformity and properties ||27|| . 

To study the properties of these media we evaluate 
bounds on the effective conductivity and elastic moduli. 
The key microstructure parameters (Ci & rji ) which occur 
in the derivation of the bounds |lj] are tabulated along 
with illustrations of the model morphologies. In addition 
we use a finite difference scheme to directly simulate the 
effective conductivity. This allows us to comment on the 
applicability of the bounds, and on their use for predic- 
tive purposes. 

The paper is organized as follows. In Section || we 



derive the 3-point correlation function for the 2-level cut 
Gaussian random field. In Section III analogous results 
are derived for the 'identical overlapping spherical annuli' 
(lOSA) or 'hollow sphere' model. In Sections |^ and ^ 
the microstructure parameters are calculated, and com- 
puter simulations of the effective conductivity are com- 
pared with the resultant bounds. In section VI we dis- 
cuss the influence of morphology on the transport and 
mechanical properties of composites. 



II. THE 2-LEVEL CUT GRF MODEL 

As in lis] we take y{r) as an isotropic Gaussian ran- 
dom field with a given field-field correlation function 
{y{ri)y{rj)) = gKinj)- Here r^ = jr^ - r^l and for con- 
venience we denote g^ [rij ) by Oii , or simply g if no am- 
biguity arises. Following Berk jlj] it is possible to define 
a composite medium with phase 1 the region in space 
where a < y(r) < /3. The remaining region is phase 2. 
In the limit /3 ^ oo the 1-level cut GRF considered in 
Rcfs. fl^,^ is recovered. The n-point correlation func- 
tion is given by 

p„(ri, r2, . . . , r„) = /f[[H{y, - a) - H{y, - /?)] \ , (1) 

where H{y) is the Heaviside function and yi — y(ji). 

The microstructure of the material is fully determined 
by specifying a, (3 and gKij"). The latter quantity is 
related to the spectral density of the field Pic(fc) by a 
clipped Fourier transform: 



gK{r) = / AT:k'^pK{k)—- — dk. 
Jo 'S'' 



(2) 



It was shown in |18| that few differences arise amongst 
the conductivity of the 1-level cut Gaussian random fields 
defined with differing spectra. Therefore we employ two 
model materials which showed the greatest variation in 
properties. In the notation of |1^ these are Model I: 

p(k) = p-V-2 ((1 _ J.2 ^ fc2)2 ^ 4^2^-1 (3) 

hm gK{r) = e , (4) 

K^oo vr 

where P is a normalization constant chosen to ensure 
3^(0) = 1, and Model III: 



pik) 



[H{p)-H{l)] (m>1) (5) 

(6) 



47r(^3 „ I 

3 (sin pr — pr cos pr — sin r + r cos r) 



r^p3 - 1) 



No normalization constant is necessary in this model pro- 
vided that K > p. In following sections we employ spec- 
trum I {ly = 0, K = cx)), spectrum I {ly = 0, K = 8) 
and spectrum III (p = 1.5). In this paper the physical 



parameters v and /x are not varied and will no longer be 
explicitly stated. 

In the notation of Appendix H the 1-point correlation 
function (or volume fraction) is just 



p = Ai(a)-Ai(/3) 



/27r 



e"2* dt. 



(7) 



The 2-point correlation function for the 2- level cut Gaus- 
sian random field can be defined as 

P2 = / + As (5, a, a) - 2A2(5, a, /3) + AaCg, /3, /3) (8) 

where we have used the fact that A2 (g, a, /3) = A2 (g, /?, a) 



and P2(.g = 0) = p . Now using Eqn. (A6) leads to 



V2{g) =/ + 



1 

2^ 



di 



- exp 



VT^^ 

1 a2 - 2a/3t + (? 

2 0^7^^ 



exp 



exp 



l + < 
1 + t 



(9) 



Similarly the 3-point correlation function is, 



P3 = P^ + A3(g, a, a, a) - A3(g, a, a, /3) - A3(g, a, /3, a) 



- A3(g, ^, a,a) 
+ A3(g,^,^,a) 



-A3(g,a,^,^)+A3(g,/3,a,/3) 
A3(g, /?,/?,/?) (10) 



where A3 is given in Eq. A13| and g = (.912, <7i3,523)- We 
could find no symmetries in these terms to allow analyt- 
ical or computational simplification of the results. 

For our purposes it is necessary to choose a and /? for 
a given value of the volume fraction -p. There are many 
ways that this can be done. An obvious method is to re- 
quire that an equivalent fraction of phase 1 lies on either 
side of a particular level cut y(v) = 7. We classify these 
'symmetric' models by the parameter 



dt 



(11) 



so that s e [0, 1]. For a given volume fraction p, a and /? 
are defined through the relations 



2 



1 



'dt 



1 



/27r 



' dt 



(12) 



Materials defined in this manner are denoted by, for ex- 
ample, III(s=.5). This indicates that the spectrum of 
model III is employed and that s — 0.5 (corresponding 
to the case a — — /3). For comparison with the 1-level 
cut case discussed in lia] it is also useful to define a 2- 
level cut GRF which reduces to the former model in a 
particular limit. This is done by fixing /? and varying 
a such that a given volume fraction is achieved. These 
'base'-level models are specified by the value 



h = 



1 



dt 



(13) 




FIG. 3. The interface of the media III(s=.5) at a vol- 
ume fraction of p = 0.2. The dark region is given by 
—0.253 < j/(r) < 0.253. Note the highly connected struc- 
tures. 




FIG. 4. The interface of the media I(s=.5) with a pore vol- 
ume fraction of p = 0.2. The light region corresponds to 
-0.253 < 2/(r) < 0.253. 



where b e [0, 1]. Since (3 is fixed, a is calculated using 
Eqn. (0). In terms of nomenclature used to describe the 
spectra previously these models are denoted, for exam- 
ple, as III(b=.3) or l(b=0) (ie. (3 = 00). The latter case 
corresponds to the 1-level cut field. 

Depending on the spectra employed and the choice of 
a and (3 the 2-cut GRF scheme can model a wide range 
of morphologies observed in physical composites. The 
morphology of 1-cut fields is characterized by a random 
array of irregular inclusions interconnected by narrower 
necks llS) similar to a 'node/bond' geometry (see 




FIG. 5. The interface of the media in(s—.2) at a vol- 
ume fraction of p = 0.2. The dark region is given by 
-1.28 < y{r) < -0.253. 



Fig. |l|(c)). This type of pore/inclusion shape has been 
observed in a range of materials including alloys Eq] and 
sedimentary rocks [£1|. Taking f3=—a in the 2-cut model 
(s = 0.5) leads to 'sheet-like' structures (see Figs. || & ^) 
with differing degrees of roughness. The smooth-sheet 
like structures of model III(s=.5) (Fig. 0) are similar 
to the pores observed in dolomitic limestone [ p9|, and 
the connected matrix in solid foams pO] (see Fig. H) and 
The rough 'sheet-like' morphology 



polymer blends |31|. 

evident in model I(s=.5) {K — 8) (Fig. ^) is similar to 
the rough porous structures observed in pore-cast stud- 
ies of sandstones |l4| , pi| . Note that certain classes of 
sandstone have been shown to have a fractal pore sur- 
face with Ds ~ 2.5 |32]. This can be modelled by taking 
iiT ^ cx) in spectrum I (Appx. H). Qualitatively different 
microstructures can be obtained in the 2-cut scheme if 
(3 y^ a. For example, the morphology of Model III(s=.2) 
(Fig. H) has both a node/bond and sheet-like quality. 



III. OVERLAPPING HOLLOW SPHERES 

A second low porosity model can be defined by gen- 
eralizing the 'Identical Overlapping Sphere' (lOS) model 
to the case of overlapping annuli (lOSA). For this model 
the probability that n points r^ chosen at random will 
fall in the void phase (ie. outside the hollow spheres) is 
just 



exp[ 



-P<")(ri,. 



0]- 



(14) 



,(") 



Here Va^ is the union volume of n spherical annuli with 

centers at r^, and p is the number density of the annuli. 

To see this consider a large region of the composite 

material of volume V which contains N = pV randomly 



positioned (i.e. uncorrelated) spherical annuli. Now con- 
sider Qn defined above. If, and only if, the center of an 
annulus is located within the volume Va\j , then one (or 
more) of the n points will lie in the solid phase. Since each 
annulus is uniformly distributed the probability that its 
center will not fall in the volume Vajj is (1 — Va^^ /V). 
Now there are iV such uncorrelated spheres so 



qn 



= lim 



pVa^ 

N 



N 



exp I 



(-p<")) (15) 



where V, and hence iV, has been taken to be infinitely 
large. This argument (for the spherical case) is due to 
W. F. Brown ||33| , p4| . By definition Qn is just the n-point 
void-void correlation function. To distinguish the cor- 
relation functions associated with the void and solid we 
refer to above model as the inverse lOSA model (as the 
correlation function corresponds to the phase outside the 
annuli). The correlation functions for the lOSA model 
(pn) are then just linear combinations of g„, Qn-i etc. 
For example, pi = I — qi and P2{ri2) = 1 — 2q + g2(»*i2)- 
Suppose the inner and outer radii of the annuli are 
p and v, then the union volume of a single annulus is 

Va = Va\j — Airiv^ — p?)/^. The number density of 
the annuli is related to the volume fraction of void (q) 
by the formula p = —\ogq/Va. The higher order union 
volumes are derived in terms of the intersection volumes 
of spheres of different radii. For the union volume of two 
annuli a distance d apart we have 



Va^^\d) = 2Va^Vaf\d) 



(16) 



(2\ 

where Vaj {d) is the intersection volume of 2 annuli. 
This function is given by 



.(2) 



A2), 



A2) 



(2), 



VaY'id) = V^tlid) + Vi^lid) - 'iVitlid) (17) 

(2'] 

with V^'{d) the intersection volume of two spheres of 
radii x and y (see Appendix O). The union volume of 
the three annuli distances a, b Sz c apart is 

Va[^^^(a, b, c) = 3Va - Vaf\a) ~ Vaf\b) - Vaf\c) 



+Vaf\a,b,c), 



(18) 



(3)n 



where the intersection volume of three annuli {Vaj ) is 



Va?'^ = V, 



(3) 
(3) 



-K 



(3) 



I^^fl 



-V, 



(3) 



+v, 



I^^v 



V, 



(3) 



I fll/fl 



luni 



^v}±.^v, 



lu^v 



(3) 



(19) 



Here the function V^J {a, b, c) is the intersection volume 
of three spheres of radii x, y and z with a the distance be- 
tween the spheres of radii y &i z, b the distance between 
the spheres of radii x &z z and c the distance between the 
spheres x & y (see Appendix y) . 




FIG. 6. The interface of the lOSA model (ro 
r\ = l.Oj at a vohime fraction of p = 0.2. 



0.8, 



As in the 2-level cut GRF model there are two obvi- 
ous ways of choosing the internal (rg = ^) and external 
(ri = v) radii for a given volume fraction. In the first 
the internal radii of the spheres is held fixed and the 
number density of spheres is increased to achieve a given 
volume fraction. This model morphology corresponds to 
manufactured materials comprised of sintered similarly 
sized hollow spheres |M . A plot of the interface for the 
lOSA model is given for the case ro = 0.8 and p = 0.2 
in Fig. pi. Using results p],p5[-p7t developed for overlap- 
ping solid spheres (ie. IOS)it is possible to incorporate 
a distribution of sphere sizes in the hollow sphere model. 
However, polydispersity effects have been shown to be 
quite small |38[| . In the second model the number den- 
sity of spheres is held fixed (so that the maximum volume 
fraction achievable is Pmax) and the internal radii is var- 
ied to achieve a given volume fraction. 

The percolation thresholds of each phase of the lOSA 
model p^ (solid) & q^ (void), can be easily derived from 
a knowledge of the threshold values of the standard lOS 
model: p^ « 0.3 |9| and g^ « 0.03 @. For the variable 
density model (ro fixed) the percolation thresholds are 
p« = 1 - (1 _p^)i-('-o/r-i)' and q^ = (<7^)i-(''«/'^i)' (so 
p^ ^ & gj^ ^ 1 as To — > ri). For the fixed density 
model the lOSA solid phase is percolative if Pmax > Pc 
and the void phase is percolative if qmin — 1 ^ 



IV. MICROSTRUCTURE PARAMETERS 

Bounds have been calculated on the conductivity [[7|JlO| 
and the bulk B and shear P,|ll| moduli of composite ma- 
terials (reviewed in Ref. ||l[J! These can be expressed |1J] 
in terms of the volume fractions and properties of each 
of the phases and two microstructure parameters: 



Ci 



2pq 



^dr 



^ds 



duP2{u)f{r,s,t) 



(20) 



5Ci 150 f°°dr f°°ds . ,„,,„, , ,„ , 

.r = ^ + ^y^ -y^ -/ duP,iu)fir^s,t) (21) 



where f{r,s,t) = pz{r, s,t) - p2{r)p2{s) /p, t^ = r^-l-s^- 
2rsu and P„ (m) denotes the Legendre polynomial of or- 
der n. As we argued in Ref. [|l8| it only appears necessary 
to know broad microstructural information about a gen- 
eral composite to successfully apply the bounds. This 
conclusion arose from the observation that the bounds 
are relatively insensitive to small variations in the mi- 
crostructure parameters. Furthermore we found that the 
parameters Ci and ryi are insensitive to fine microstruc- 
tural details within a class of composites (e.g. the over- 
lapping sphere class, or the 1-level cut GRF class). An 
example of this insensitivity is also seen when polydisper- 
sity effects of particulate models are considered pS] . In 
light of these facts the parameters calculated from mod- 
els may well have application to physical composites for 
which precise microstructural information is unavailable. 

In Fig. |7| we provide a graphical summary of the wide 
range of isotropic composites for which p^ (and hence 
the microstructure parameters) can been exactly calcu- 
lated. It is clear that the 2-cut GRF and overlapping 
hollow sphere model considerably expand the classes of 
materials to which the bounds can be applied. 

We now report calculations of the microstructure pa- 
rameters for a variety of 2-level cut GRF and lOSA mod- 
els. Our method of calculating C,i (and 771) has been dis- 
cussed previously |1^ . In addition we employ an adaptive 
integration algorithm |^1|] to compensate for the fact that 
the sub-integrand J_^Pn{u)f{r,s,t)du varies rapidly in 
the region r « and involves a considerable number of 
function evaluations. The error in the results is less than 
1%. To model as wide a range of materials as possible 
three qualitatively different spectra are used in the level- 
cut GRF scheme; models I {K = 00), I {K = 8) and III 
{fi = 1.5). These spectra lead to surface-fractal, rough 
and smooth interfaces respectively. 

As we are primarily interested in low volume fraction 
porous or solid media the microstructure parameters we 
report are in the range 0.0 < p < 0.4. The results for (i 
and ?7i are given in Tables 13 & and selected results are 
plotted in Figs. || & y. The results for the two variants of 
the lOSA model are given in Table pl| and plotted along 
side the results for the 2-cut GRF models in Figs. || & ||. 
Due to the simple geometry of the lOSA model it is pos- 
sible to calculate Ce to order p (see Appendix O). This 
result can then be used to show Ci|p=o = (^o/'^i)'^ (rep- 
resented by symbols in Fig. M) in agreement with our 
numerical calculations of Ci . 

To compare the properties of different media we plot 
(Fig. no) the upper bound on the conductivity for one 
member of each class of composite: 2-cut GRFs, hollow 
spheres, lOS-voids |33 (or swiss-cheese) , 1-cut GRFs IT^] 
and lOS 134] (or solid spheres). 
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(a) lOS, T=8 



(b) lOS voids, T=8 



(c) IPSA, ro=0.8. T=8 
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(f)I(b=0),v=4,K=16, T=47l 




(g) I(s=.2), T=47l 






(h) I(s=.35), T=47l 







(j) III(s=.2), T=87t (k) in(s=.35), T=87t (I) in(s=.5), T=87t 

FIG. 7. A variety of microstructures (volume fraction 20%) 
for which the third order statistics are known exactly: (a)-(c) 
lOS and lOSA models (n = 1); (d)-(f) 1-level cut GRF's; 
(g)-(l) 2-level cut GRF's. Except where noted Model I has 
K = 8 &u = and Model III /j, = 1.5. 



TABLE I. The microstructure parameter C^i for a range of 
materials generated from the symmetric GRF modcL 



Mod. 


I, 


K = oo 


I, 


K = 


B 




III 




s 


.20 


.35 


.50 


.20 


.35 


.50 


.20 


.35 


.50 


P 










Ci 










.050 


.401 


.401 


.402 


.706 


.773 


.786 


.785 


.872 


892 


.075 


.402 


.408 


.409 


.641 


.719 


.739 


.733 


.845 


873 


.100 


.405 


.413 


.415 


.597 


.684 


.706 


.691 


.824 


858 


.125 


.410 


.422 


.425 


.563 


.655 


.677 


.655 


.807 


845 


.150 


.417 


.428 


.431 


.536 


.633 


.656 


.625 


.791 


828 


.200 


.425 


.443 


.449 


.500 


.601 


.628 


.574 


.769 


819 


.250 


.435 


.459 


.464 


.478 


.583 


.611 


.532 


.753 


811 


.300 


.443 


.475 


.481 


.464 


.575 


.605 


.495 


.741 


808 


.350 


.451 


.491 


.497 


.455 


.572 


.603 


.456 


.734 


807 


.400 


.456 


.506 


.515 


.444 


.574 


.607 


.411 


.728 


810 



TABLE II. The elasticity microstructure parameter ryi for 
a range of materials generated using the symmetric GRF 
model. 



Mod. 


I, K = oo 


I, K = 8 


III 




s 


.20 .50 


.20 .50 


.20 


.50 


P 




Vi 







050 


.355 


.351 


.523 


.613 


.609 


.754 


075 


.358 


.362 


.463 


.548 


.543 


.705 


100 


.362 


.369 


.430 


.516 


.500 


.672 


125 


.370 


.377 


.416 


.493 


.471 


.648 


150 


.373 


.388 


.407 


.480 


.449 


.608 


200 


.394 


.402 


.404 


.473 


.426 


.621 


250 


.410 


.430 


.410 


.478 


.414 


.609 


300 


.426 


.451 


.420 


.492 


.408 


.615 


350 


.438 


.474 


.430 


.510 


.406 


.627 


400 


.442 


.495 


.431 


.533 


.396 


.643 



TABLE III. The microstructure parameters for different 
versions of the lOSA model. For the case ro = the results 
are just those of the standard lOS model (see Ref. ||lj) 



Mod. 


''o = 


= 0.5 


ro = 


= 0.9 


Pmax 


= 0.9 


Pmax 


= 0.7 


P 


Ci 


Vi 


Ci 


Vi 


Ci 


Vi 


Ci 


Vi 



05 


.152 


.119 


.737 


.468 


.974 


.936 


.955 


.888 


10 


.179 


.153 


.744 


.490 


.948 


.880 


.911 


.788 


15 


.207 


.187 


.752 


.512 


.924 


.827 


.870 


.707 


20 


.234 


.221 


.759 


.533 


.900 


.780 


.829 


.640 


25 


.262 


.254 


.766 


.554 


.877 


.746 


.791 


.590 


30 


.289 


.288 


.772 


.576 


.856 


.710 


.754 


.551 


35 


.317 


.322 


.778 


.596 


.835 


.683 


.718 


.524 


40 


.345 


.356 


.784 


.616 


.815 


.662 


.683 


.508 


50 


.402 


.424 


.794 


.656 


.776 


.638 


.614 


.502 


60 


.459 


.492 


.801 


.697 


.741 


.631 


.539 


.520 


70 


.517 


.560 


.805 


.733 


.706 


.643 


.414 


.511 


80 


.578 


.630 


.804 


.771 


.666 


.668 






90 


.643 


.705 


.792 


.804 


.558 


.658 







We have also evaluated bounds on the shear, bulk 
and Young's moduli of the models. In Ref. |30| we 
showed that the upper bound on Young's modulus was 
in good agreement with experimental measurements for 
foamed solids. Model III(s=.5) provides a good model of 
polystyrene foam (compare Figs. || & ^, and the lOSA 
model accurately mimics the microstructure and proper- 
ties of sintered hollow glass spheres. In Fig. |lj the upper 
bound on the shear modulus is shown for each class of 
composite considered above: the microstructure clearly 
has a strong influence on elastic properties. The bulk 
and Young's moduli show similar behaviour. 

V. SIMULATIONS OF ae 

In addition to bounding the properties of composite 
media and providing qualitative information on these 
properties, it has been observed that the bounds also 
have reasonable predictive power ffl . To test the predic- 
tive utility of the bounds and provide a direct 
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FIG. 8. The microstructure parameter ^i for selected mod- 
els. The lOS model and model III(h=0) are included to show 
the behaviour ofC^i for different classes of composites (see Fig. 
M). The solid symbols represent analytic calculations of ("i|p=o 
for the lOSA model. 
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FIG. 9. The elasticity microstructure parameter rji for se- 
lected models. The lOS model and model III(b=0) are in- 
cluded to show the behaviour of rji for different classes of 
composites (see Fig. m). 

comparison between microstructure and properties we 
use a finite-difference method to explicitly calculate the 
conductivity of several 2-cut GRF's. 

The effective conductivity a,, of a composite is defined 
as the ratio of the current density to the applied poten- 
tial. We take T as the scale of the sample and M^ as 
the number of nodes (so the spatial resolution scale is 
Aa; = T/M). The generation of random fields and the 
method for determining ae were described in Ref. [|l8| for 
the case of 1-cut fields. A number of additional difficulties 
are encountered in the simulations of Ue for the 2-level 
cut GRF's. The major problems are; (i) discretisation 
effects which occur when the discretisation length scale 
Aa; is insufficient to resolve the thin sheet-like structures 
which arise (e.g. Fig. 0) and (ii) finite-scale effects which 
arise if T is not large enough to represent an 'infinite' 
medium. In practice T should be several times the 
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FIG. 10. The 3-point upper bound M on the effective con- 
ductivity (contrast 1:0) of five different classes of microstruc- 
ture. The data for the lOS and 1-cut GRF models are from 
Refs. J3i|ii 
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FIG. 11. The 3-point upper bound jJlf on the effective 
shear modulus (contrast 1:0) of five different classes of mi- 
crostructure. The data for the lOS and 1-cut GRF models 
are from Refs. Edili 



correlation length of the microstructure (approx. unity). 
Discretisation effects can be reduced by increasing M or 
decreasing T (to increase the width of the sheets relative 
to Ax). However our computational requirements dic- 
tate M < 128 and decreasing T leads to noisy results. 
Thus T must be chosen to minimize each of these com- 
peting errors. By performing several numerical tests |4l| ] 
a reasonable value of T was determined to ensure that 
simulations of Cp are accurate. As the sheets become 
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FIG. 12. The simulation data and bounds for three differ- 
ent 2-level cut GRF models at contrast cri,2 = 10, 1. The 
thick solid line corresponds to the bounds for the 1-level cut 
model III(b=0). The bounds are clearly seen to differentiate 
the different classes of media at low p. 



FIG. 13. The simulation data and bounds for three differ- 
ent 2-level cut GRF models with conductivity ai,2 = 50, 1. 

TABLE V. The effective conductivity of several 2-level 
cut GRF's for the case ai 2 = 50, 1. 



TABLE IV. The effective 


conductivity 


of several 2-level 






III(s=. 


2) 


III(s=.5 


) 




I(s=.5) 




5 lor Lue 


case 


cri,2 - 


- 10, 1. 










P 
0.05 
0.10 


T/tt 
4 
4 


0"e 

2.3 
3.4 


Err. 
0.2 
0.4 


T/n 
2 
2 


2.6 
4.1 


Err. 
0.2 
0.4 


T/tt 
1 
1 


o"e Err. 






III(s=.2) 




III(s=.5) 






I(s=.5) 




2.1 0.1 


P 


T/ti 


ffe 


Err. 


T/tt 


(Te 


Err. 


T/n 


O-e 


Err. 


3.3 0.1 


0.05 


4 


1.28 


.01 


2 


1.32 


.02 


2 


1.27 


.00 


0.15 


8 


4.1 


0.1 


2 


5.7 


0.5 


1 


4.2 0.1 


0.10 


4 


1.51 


.02 


2 


1.64 


.03 


2 


1.53 


.02 


0.20 


8 


5.1 


0.2 


2 


7.4 


0.7 


1 


5.4 0.1 


0.15 


4 


1.76 


.04 


2 


1.96 


.05 


2 


1.77 


.02 


0.25 


8 


6.4 


0.2 


4 


9.3 


0.2 


2 


7.0 0.1 


0.20 


4 


2.02 


.06 


4 


2.30 


.07 


2 


2.05 


.01 


0.30 


8 


7.8 


0.2 


4 


11.2 


0.3 


2 


8.6 0.1 


0.25 


8 


2.32 


.03 


4 


2.58 


.10 


4 


2.35 


.01 


0.35 


8 


9.3 


0.3 


4 


13.2 


0.3 


2 


10.4 0.2 


0.30 


8 
8 


2.61 
2.93 


.04 
.03 


4 
4 


2.94 
3.30 


.12 
.13 


4 
4 


2.67 
3.02 


.01 
.01 


0.40 


8 


10.9 


0.4 


4 


15.3 


0.3 


2 


12.3 0.2 


U.35 




















0.40 


8 


3.25 


.05 


4 


3.70 


.13 


4 


3.40 


.02 





















thinner (ie. p decreases) it was found that smaher val- 
ues of T are necessary to ehminate finite-scale effects. 
This can be explained in terms of the faster decay of the 
correlations between the components of phase 1. 

We choose to study the effective conductivity of mod- 
els III(s=.5), I(s=.5) (K = 8) and III(s=.2). The former 
models provide examples of smooth and rough 'sheet- 
like' pores. The latter model (III(s=.2)) has a mor- 
phology comprised of inclusions with both a sheet- and 
node/bond-like quality. The conductivity contrasts em- 
ployed occur in physical composites and have been stud- 
ied previously, allowing comparisons to be made. In each 
of the cases we report results averaged over five samples 
for a range of volume fractions 0.05 < p < 0.4. In all 
cases the simulational data lie between the bounds. 

First we consider the conductivity contrast ai^2 — 
10, 1. The results are tabulated in Table [V, and plotted 
in Fig. n2 along with the bounds for each model, Here, 



and in subsequent calculations, the lower bound of Be- 
ran and the upper bound of Milton (see Ref. [|8|) 
are employed. The data for model III(s=.5) practically 
lies along the relevant upper bound. In contrast the ef- 
fective conductivity of models III(s=.2) and I(s=.5) fall 
between the bounds, however the upper bound still pro- 
vides a reasonable estimate of CTc in each case. For pur- 
poses of comparison the bounds for the 1-level cut GRF 
model III(b=0) are included in Fig. n2. At low p the 
bounds clearly differentiate between the different classes 
of media. It is clear that the model III(s=.5) is a signifi- 
cantly more efficient conductor than models III(s=.2}or 
I(s=.5) and those defined using a 1-level cut GRF in [[18[ . 
The simulation data for the contrast (Ti_2 = 50, 1 is re- 
ported in Table ^ and plotted in Fig. |lj. Qualitatively 
the results are the same as those discussed in relation to 
the case ai^2 = 10, 1. Note that the upper bound is again 
a good estimate for model III(s=.5); less so for models 
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FIG. 14. Simulation data and bounds for conductivity con- 
trast CTi,2 = Ij 0. The models and contrast are relevant to 
solid-foams /pOl/ and porous rocks MTJI. 

TABLE VI. The effective conductivity of several 2-level 
cut GRF's for the case ai? = 1, 0. 



FIG. 15. Simulations of a^ and the bounds for cri,2 = 0, 1. 
This contrast is relevant to modeling transport in membranes. 



TABLE VII. The effective conductivity of several 2-level 
cut GRF's for the case ai 2 ~ 0, 1. 
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III(s=.2) and I(s=.5). 

In porous rocks and solid foams the conductivity of the 
medium surrounding the conducting pathways has negh- 
gible (or zero) conductivity. To model such systems the 
contrast cri^2 = 1, is used. The data and computational 
parameters used in the simulations are reported in Table 
VI. Each material is seen to be conductive at the lowest 
volume fraction considered p = 0.05. Discretisation ef- 
fects prohibit accurate simulations of (7^ at lower volume 
fractions. The simulation data and the upper bounds 
are plotted in Fig. |lj. Even in this large contrast situa- 
tion the upper bound for model III(s=.5) agrees with the 
data. 

To consider the case of diffusive transport in mem- 
branes, we assume that the membrane has negligible dif- 
fusivity with respect to the surrounding fluid. Therefore 
the contrast ai^2 = 0, 1 is employed. For this system large 
discretisation effects prohibit the consideration of model 
I(s=.5) and membrane/pore volume frac tions of less than 
p = 0.10. The data is presented in Table VII and Fig. UH. 



Note that the presence of a membrane occupying 10%- 
20% of the total volume reduces the diffusivity by a fac- 
tor of two. This is due to the tortuous pathways through 
which the diffusing species must migrate. In contrast to 
three cases considered above the upper bound does not 
provide a good estimate of o-g for model III(s=.5). 

VI. EFFECT OF MICROSTRUCTURE ON 
PROPERTIES 

The precise role of microstructure in determining the 
macroscopic properties of composite media has been the 
subject of many studies. A number of simple models 
of pore-shape have proposed to determine, for exam- 
ple, the effect of pore-size distribution [^, pore rough- 
ness [Q and pore geometry |2^ on transport in porous 
rocks. Simple micro-mechanical models |45(] have also 
been studied to ascertain, for exa mple , the effect of inclu- 
sion shape [Q and cell structure |Q,Q on the mechanical 
properties of composites. In this section we investigate 



TABLE VIII. Qualitative comparison of microstructure and macroscopic properties. We consider the contrast 
1:0 at a representative volume fraction p = 0.2. Uu, k„, and /!„ are upper bounds on the effective conductivity M, 
bulk M and shear modulus [hll respectively. 
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how morphology influences the properties of realistic 
model composites. 

To simplify the discussion we summarize relevant data 
for a variety of G RF and particulate microstructure mod- 
els in Table VIII, We consider systems of 1:0 contrast at 
p = 0.2. This case corresponds to a conducting (mechani- 
cally strong) matrix in an insulating (weak) medium (e.g. 
foamed solids l^^l ) ■ This contrast also corresponds to low 
porosity conducting pores in an insulating medium (e.g. 
porous rocks). To gauge the effect of microstructure on 
material properties we assume that the upper bound on 
each property provides an estimate of its actual value. A 
comparison of cr„ and cr^ (where available) shows that this 
is generally true l^^l- Note, however, that if the difference 
between au for each of the models is small (e.g. IPS- voids 
and model III(s=.2) and examples in Ref. |18]) such an 
assumption cannot be made (cr^ > c^ but ct" < a^). 

At a 1:0 contrast the effective properties only differ 
from zero if the composite is macroscopically connected 
(i.e. percolative). At p = 0.2 this condition is satis- 
fied for all but one of the media (conducting spheres in 
an insulating medium). Above this threshold the mag- 
nitude of the macroscopic properties is then governed 
by the shape of the inclusions. It is clear from the ta- 
ble that sheet-like structures provide higher conductivity 
and mechanical strength than those with a node/bond- 
like character. To elucidate the role of inclusion shape 
we derive approximate expressions for the effective con- 
ductivity of periodic media with unit cells of each type 
in Appendix H. For small volume fractions (p <C 1) the 
node/bond model has Ue ~ p^cri and the sheet-like model 
has CTe — ■jPO'i ill qualitative agreement with the data. 
Interestingly the periodic sheet model provides a surpris- 
ingly good estimate of CTe = §0.2 = 0.133 for model 
III(s=.5) (o-e = 0.130). 



From the table it also evident that interfacial roughness 
plays an important role in determining properties. Con- 
sider the 2-cut fields I(s=.5) and III(s=.5). In Figs. |(i) 
and |j(l) it is clear that both models contain sheet-like 
pores. The differences are then due to the interfacial 
roughness. This is confirmed by comparing the relative 
values of ct„ for model I in the cases K — 8 (smooth 
on scales below Xmin — 2tt/8) and K — oo (rough on 
all scales) - see Appx. ^. The effect of increasing K 
from 8 to 32 on the morphology of model I(s=.5) can be 
seen by comparing Fig. 7|(i) & Fig. |l^(a). In the rough 
model the sheet-like pores are thinner and a large pro- 
portion of pore space is distributed in protrusions. As 
these protrusions contribute little to the overall conduc- 
tance (or strength), this significantly reduces both con- 
ductivity and strength. This also explains why Model III 
is more conductive (stronger) than Model I. The much 
smaller effect of roughness on morphology of 1-level cut 
fields can be seen by comparing Fig. ^(d) & Fig. p^(b): 
the basic inclusion shape is less affected than in the 2-cut 
case. 

Now cons ider t he data for the spherical particulate me- 
dia in Table VIII . The hollow sphere model appears to be 
more conductive, or stronger, than the lOS-voids (swiss- 
cheese) model. This is due to the fact that the former 
model has an approximately sheet-like character (Fig. g) 
in contrast to the node/bond-like structures H apparent 
in the inverse lOS model. The lOS model at p = 0.2 
does not have a sufficient density of spheres to provide a 
percolative pathway. 
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VII. CONCLUSION 



APPENDIX A: THE LEVEL CUT GAUSSIAN 
RANDOM FIELD 



In this paper we have derived the 3-point statistical 
correlation functions for two models of random compos- 
ite media. The results were applied in the evaluation of 
bounds on the effective conductivity and elastic moduli 
of each model. In addition the 'exact' effective conduc- 
tivity was estimated for the 2-level cut GRF model by 
direct simulation. The models are applicable to physical 
composites which remain percolative at very low volume 
fractions Pc < !%• These include solid foams, porous 
rocks, membranes and sintered hollow glass spheres. In 
contrast, previously employed models of microstructure 
have percolation thresholds of order 10%. 

Microstructure was demonstrated to have a strong in- 
fluence on the effective properties of composites. The 
relative variations amongst the 2-level cut, lOSA and 1- 
level cut models were attributed to three morphological 
factors; pore-shape, interfacial roughness and the perco- 
lation threshold of the material. Materials with sheet-like 
inclusions were shown to have a significantly greater con- 
ductivity/strength than materials with node/bond-like 
inclusions. By comparing the microstructure parameters 
of similar composites with fractally rough and relatively 
smooth inclusions we found that interfacial roughness de- 
creased composite conductivity/strength. The observa- 
tion was confirmed by directly comparing simulated val- 
ues of CTe for Model I(s=.5) (sheet-like and rough) and 
Model III(s=.5) (sheet- like and smooth). The behaviour 
was physically attributed to the fact that the protrusions 
of rough interfaces contribute httle to effective properties. 

The models discussed here considerably expand the 
range of systems to which bounds can be applied. To 
facilitate use of these bounds we have tabulated cross- 
sections and microstructure parameters for a number of 
different variants of each model. Such bounds have two 
clear applications. Firstly they can be used to narrow the 
possible microstructures of a composite for which prop- 
erties are known; composite materials may violate the 
bounds for a particular model system. Indeed for certain 
cases of realistic media the bounds are mutually exclusive 
(see Fig. |2|). Secondly, the upper bound is often a very 
useful estimate of the actual property. Indeed for model 
III(s=.5) the upper bound provides an excellent estimate 
of the effective conductivity over the full range of volume 
fraction measured. 

The simulational data presented here allows compari- 
son of model properties with those of physical compos- 
ites [^^lii ill- Furthermore the data can be used to 
assess both predictive theories for ct^, and higher order 
bounds. We note that the 4-point correlation functions 
of the two models considered here can be calculated, and 
hence used to evaluate known 4-point bounds ||ll|,Q . Fi- 
nally we remark that the generalization of the lOS model 
to include the case of hollow spheres broadens the utility 
of the model as a 'bench-mark' theoretical tool; as well 
as providing a realistic model of certain composites. 



In this appendix several results are derived which 
are useful for calculating the n-point correlation func- 
tion of a material defined by level cut(s) of a Gaussian 
random field |l^,^,^. The joint probability density 
(JPD) of a Gaussian random field is -Pn (2/12/2 •■• ?/«) = 
((27r)"|G'|)-5 exp(-iy^G'-V)- where the elements of G 
are 5^- = g{nj) = {y{ri)y{rj)) Q. The function g has 
the properties 17(0) = 1 and limr^oo g{f) -^ 0- By defini- 
tion we have p„ = /^ dyi . . ./^ dynPniviVi ■ ■ ■ Vn)- Note 
that, in this form, p„ is difficult to evaluate. For example 
if gij ~ 1 for all i,j then \G\ ~ (n > 1). It is possi- 
ble to avoid such problems, and reduce the number of 
integrations required, by taking the following approach. 

Expanding Eqn. (nh gives terms of the form. 



A„(g,a) = / J|iJ(y, -a,) 



(Al) 



\i=i 



where g = (gi2,. .. ,5(«-i)«), a = (ai,...,a„) and the 
ai are equal to a or /?. The analysis which follows relies 
on an integral representation of the Heaviside function. 



-1 



H(y- a) = ^r— I e 



dw 



27ri ,/c ^ 



(A2) 



where the contour C lies along the real axis except near 
the origin where it crosses the imaginary axis in the upper 
half plane. 



Now we turn to the evaluation of the terms Eqn. ( Al) 
For the case n = 1 we have G = gn = 1 so 



Al = (iJ(yi - ai)) = 



1 



dt. 



(A3) 



Now consider A2, in this case the matrix G in the JPD 

is 



G 



.911 .912 




[1 .9] 


.921 .922 




[,9 ij 



(A4) 



with g — gi2 — 321 and \G\ — \J\ — g'^. Using the Heav- 
iside function and interchanging the order of integration 
gives. 



A, 



1 /■ dwx f dw2 i^T 

-e 
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X / dyi I dy2 e" 
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In this case we differentiate with respect to g 
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^=-7^p-ip-.e-^^^«^+-^-, (A5) 
dg (27ri)2y-oo J-oo 

and perform the integrals with respect to Wi. The result 
is simply integrated to give A2 (up to a constant) 



A2 = 



1 



(2^)70 VT^^ 



dt ( a? — 2aia2t + Qfo 

exp ' 



2(1 -t2) 



(A6) 



The derivation of A3 follows similar lines: The initial in- 
tegration over the yi gives 



A, 



-1 



dw\ f dw2 f dwj, 



(27ri)3 J (J wi Jc W2 Jc W3 



-e 2 



w Gw+2\v a 



(A7) 



For this case |G| = 1 - 5^3 - 9i3 - ffis + ^912913923 
and w^Gw = wf + W2 + w^ + 2wiW2gi2 + 2wiW3(?i3 + 
2w2W3g23- Taking the derivative of Eqn. (A7) with re- 
spect to gi2 gives 
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where w = (^1,102), u = (mi,U2) = ("i + igi3W3,a2 + 
ig23W3) and 
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1 512 
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Performing the standard integrals with respect to Wi and 
W2 gives, 
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2-KiJ Jc W3 

where i' = |G|/(1 — 5^2) ^'^^ 

ai(ffi3 - 912923) - "2(323 - 512513) 

K = ^3 2 ■ 

1-512 

Now the remaining integral can be re-expressed to give 



(AlO) 



(All) 



1 



1 



aA3 

dgi2 (27r) ^/fZ 

1 



9h 



exp 



aj - 2aia2gi2 + a\ 



2(1 -5?2) 



1 ^2 , 



(A12) 



i^. 



|G| 



X au 



ai[gik - 9jk9ij) + aj{9jk ~ 9ik9io) 

1-4 



The results can be formally integrated to give, up to a 
constant, 

A3(g,a)= / dt[gi2Ai2{t^,a) 
Jo 
+913 Ai3 (tg, a) + 523^23 (ig, a)] . ( A13) 

The results for A^ are employed in the text to derive the 
statistical correlation functions. 



APPENDIX B: FRACTAL SURFACE 
DIMENSION 

Berk [g5[ has shown that the class of level-cut GRF 
models with spectra p{k) ~ j^k^^'^^^ as /c -^ 00 (0 < 
e < 1) have field- field correlation functions (7 (r) ~ 1 — br^'^ 
and surface fractal dimension Dg = 3 — e. Here a and 
b are related constants. In this appendix we show how 
the finite cut-off wave-number K effects the roughness 
(fractal) properties of a GRF interface. Through a very 
elegant argument Debye et al [Q showed that the sur- 
face to volume ratio {S/V) of a porous solid was related 
to the two point correlation function by 



V2(0) 



s_ 

V' 



(Bl) 



Now consider p2 for the general 2-level cut Gaussian ran- 
dom field. The most instructive method of examining 
^2(0) is by generating an expansion for small r. Thus we 
write 



P2{r) 



1 



dt 



^^ 27r J,_s^T~P 



fit), 



(B2) 



where d{r) = 1 — g{r), and f{t) is a suitably defined 
function. Integrating by parts and retaining leading or- 
der terms gives 



P2{r)^p-y/2S(^)f{l)/27r, 



(B3) 



with /(I) = exp(-aV2) + exp(-/3V2). Now if S{r) = 
0{r^) then the specific surface is well defined and ^2(0) 
can be evaluated. However for the class of spectra con- 
sidered by Berk |25| 6{r) ~ br^'^ so 



where F12 = nj ^. Similar expressions can be de- 
rived for 9A3/9512, and dh.3ldg23- These are denoted 
by Aij(g,a) = dh.3/dgij. With k j^ i 01 j we can also 
write a general expression for Fij 



P2^P^ 



26/(l)r727r. 



(B4) 



Therefore ^2(0) ^^'^ the specific surface (S/V) are infi- 
nite. Bale and Schmidt [^ have shown that this type 
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(a) I(s=.5), K=32, T=4jc 



(b) l(b=0), K=32, T=47C 



FIG. 16. Cross sections of two models based on spectrum 
I. This figure shows the roughness of the interface for large K 
(Compare Figs, ^(i) & (d)). 

of singular behaviour implies a fractal surface. The frac- 
tal dimension 'Z3s' is given in terms of the correlation 



function through the relation p2 (r) ~ p ^ 



„3-D^ 



with 



some constant. We infer from Eqn. (B4) that our Model 
I (e = 1/2) has a fractal surface with Dg = 2.5. 

As discussed in W^ it is necessary to introduce a finite- 
cutoff wave number K for computational and physical 
reasons. We now show how this parameter changes the 
microstructure. The wave-number K corresponds to a 
cut-off wavelength Xmin = '2tt/K which specifies the scale 
of the smallest "ripples" on the surface. Thus we expect 
the surface area to scale as a fractal down to some length 
scale related to Xmin- This can be confirmed mathemat- 
ically by considering the small r behaviour of p2 (5k (?"))• 

For arbitrary r, QKi^) can be expressed in terms of the 
moments of k. Using a Taylor series expansion for sinfcr 
in the definition of g^ (0) we have 



9K{r) 




47rp,^(fc)r<ifc r^ + O(X'V^) 



where the latter approximation is valid if r ^ X mm ■ Sub- 
stituting this result into the expansion for p2 (B3) and 
using relation (Bl) gives p^,p6[ 



5 
V 



!#(-^' 



-ip' 



(B5) 



Thus for r <C Xmin the surface is behaving in a regular 
manner {Dg = 2) as anticipated. Note that for the case 
K —> oo and e < 1 the moment (fc^) diverges and this 
approximation does not apply. 

To examine the behaviour for r > Xmin we can succes- 
sively integrate (0) by parts |18| to obtain 



9 Kir) 



9oo{r) 



(B6) 



pX2e \ ^2 



COS Kr 3 sin Kr 



^3^3 



0{K-S 



If Kr > 1 this expansion is asymptotic to g^ IM . Now 
in the region Xmin <C r ^ 1 the algebraic terms in the 



expansion are negligible and g^ — g — {l — br'^'^)/P (with 
P^l). 

In summary we have 



P2(r) 



i(fc2)m))^ 0<r«A„ 



/(I) 

27r 



r' Xr 



< r< 1. 



(B7) 



This demonstrates the regular [Dg — 2) nature of the 
surface in the former region, and the fractal behaviour 
{Ds ~ i — €) over the spatial scales in the latter region. 



APPENDIX C: INTERSECTION VOLUME OF 
TWO AND THREE SPHERES 



The intersection volume V^ ' id) of two spheres of radii 
^ and V separated by a distance d is simple to calculate. 
With ri = min(/i, i^) and r2 — max(/i, i^) Vj — 47rrJ/3 
if < rf < r2 - ri , y/^-' = if r2 + ri < d < oo and 



Vl^'^ = 2nirl + r3)/3 - n{rjxi + rlx^ 



A) 



(Cl) 

if r2 - ri < d < r2 + ri. Here xi = {d^ + rf - r2)/(2d) 
and X2 = d — xi. 

A compact form of the intersection volume of three 
spheres of equal radii [r — 1) has been derived previ- 
ously by Powell [^ . Several of the key simplifications in 
the derivation formula are not possible when the spheres 
have different radii. However a less elegant but straight 
forward result can be determined. Suppose the spheres 
have radii va, tb and re and are distances a, h and c 
apart and that there exist two unique points P and Q 
where the surface of the spheres meet. From Powell M] 
the intersection volume of the three spheres is equal to 
twice the following expression (Powell's Theorem): 

The volume of the tetrahedron PABC 

— The volume of the sphere center A enclosed by the 
faces of the tetrahedron PABC which meet at A 

— The volume of the sphere center B enclosed by the 
faces of the tetrahedron PABC which meet at B 

— The volume of the sphere center C enclosed by the 
faces of the tetrahedron PABC which meet at C 

+ The intersection volume of the spheres centered at B 

and C enclosed by the two faces of the tetrahedron 

PABC which meet in BC 
+ The intersection volume of the spheres centered at C 

and A enclosed by the two faces of the tetrahedron 

PABC which meet in CA 
-|- The intersection volume of the spheres centered at A 

and B enclosed by the two faces of the tetrahedron 

PABC which meet in AB 
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(0,0) ^ 




aB 


(a,( 


c\^ 






/ 


bc^\ 


(xp,yp) , 




/S 



A (b cos C,-b sin C) 



(b) 




FIG. 17. Geometrical elements used in the calculation of 
the intersection volume of three spheres of different radii. The 
axis and nomenclature are those employed in the text, (a) 
the tetrahedron base; (b) the sphere centered at A enclosed 
by three faces of the tetrahedron. 

The cases where there is no unique point of intersection 
between the spheres is discussed below. We first define a 
coordinate system with origin at the center of sphere C 
as drawn in Fig. |l^(a). By solving the equations of the 
three spheres simultaneously it is simple to show that 



Xp = 



yp = 



rs 



re 



Zp 



-b^ + 


2a 

„2 


-rl + 2b 


COS 


Cxp 


\Ac- 


Xp 


2bsmC 
Vp- 







(C2) 
(C3) 
(C4) 



It is also necessary to know the distances aB,ac ■ ■ ■ given 
in Fig. |l^(a). We have ac = (a^ + ^c ^ ^|)/(2a), 



rl)/i2b) CB 
- b — be and ca 



K +r%- rl)/(2c), 



CB- 



bA = {V 

aB = a- ac, bA 

The volume of the tetrahedron is W = ^absinCzp. 
The solid angle (J)a of the tetragonal wedge at A (see 
Fig. |3(^)) ^^^ ^^ calculated by using the fact that 
(l)A = (E + F + D - n) and 



cosD 



cos d — cos e cos / 
sin e sin / 



(C5) 



(similarly for cos E and cos F) . This gives 




FIG. 



(c) (d) 

18. Several of the topologically distinct cases which 
arise in the calculation of the intersection volume of three 
spheres of different radii. 



COS 



Ca — cos Ab/ 



smAy/r\ — 53 



cos A — bACA 



x/r\-b\y^ 
b^ — cos Ac A 



sinyl^r^- 



(C6) 
(C7) 



Similar results are obtained for the solid angles 4>b.c- 

It is critical to know whether the point xp,yp lies in- 
side or outside each of the faces of the triangle. This can 
be done by defining the variables 



SA = sgn{-yp) 

SB ~ sgn(cos Cyp + sin Cxp) 

sc — sgn(cos Byp — sin Bxp 



asinB) 



(C8) 

(C9) 

(CIO) 



Then for example sa = ±1 as the point {xp,yp) is in- 
side or outside face a of the triangle ABC. In the case 
rA = rB = re ~ I (Powell ||5|]) we have xp — a/2 
and yp — — ccos^/2sinC so that sa — sgn(cos^), 
SB = sgn(cos B) and sc — sgn(cos C) as they should. 

The wedge angle associated with the intersection vol- 
ume of spheres B & C is, 



0A 



SA 



Vr% -a%- 



(Cll) 



Similarly for the angles 9b ^ dc- 

Now the volume of a tetragonal wedge of solid angle 
is r^(f)/'i and the intersection volume of spheres enclosed 



14 



in a wedge of angle 9 is 9Vj^ ' /2n. Therefore, by Powell's 
theorem, 



VifyM^b,c) 



2Vt 



r 



9A 



3^ 



9B 



^z'ct>c (C12) 



IT " TT IT " 



Here x = ta, y — tb &^^ z = re- This formula is equiva- 
lent to Powell's result |53| in the case x — y ~ z — 1. 

Several other cases arise if the point P does not exist. 
Some of these are illustrated in Fig. UH. Either two (or 
more) of the spheres are disconnected (not illustrated), 
they are connected but Va\ = (b) or the intersection 
volume is given by that of two of the spheres (c) or some 
other formula (d). 

APPENDIX D: DERIVATION OF Ci p=o (IOSA) 

It is possible to develop an independent check on the 
calculation of (i for the IOSA model by direct calculation 
of CTe ■ Using the framework of Reynolds and Hough ]5^ ] 
gives 



Ce = 0-2 + (o-i - o-2)p/i 



(Dl) 



where /i — Ei/E. Here Ei is the average of the field 
throughout phase 1 and E is the applied field. While the 
above formula is exact it is only possible to evaluate Ei 
approximately. In the low concentration regime (p ^ 1) 
El is the field within a hollow sphere (conductivity cri) 
embedded in an infinite medium (conductivity (72) sub- 
ject to an applied field E. To determine this field we con- 
sider a more general problem where the conductivities of 
the innermost spherical region (0 < r < tq), the annulus 
[tq < r < ri) and the enclosing medium (ri < r < cxi) 
are Oa-^ o^ and a^ respectively. The potential of the 
field satisfies Laplace's equation and charge conservation 
boundary conditions at phase boundaries. Using stan- 
dard techniques it is possible to show that, in each region, 
the potential has the form (f)d — {Adv + Bdr^'^) cos 6* with 
d = a,b OT c. Applying the appropriate boundary condi- 
tions on each of the faces of the hollow sphere gives. 



Aa 

Bb 
Be 



^9s^E/H- Ba^O; Ab 
-3rlE{l-x)/H- Ac^-E; 

-rlE[{l-x){l + 2y) + s^2- 



3s^E{2 + x)/H; 
x){l-y)]/H. 



where H = 2{1- x){l-y) + s^{2 + x){2 + y), x = cja/ab, 
y — cTb/ac and s — rijrQ. For the desired value of /i, 
Ca = Tc = (72 and Ob — oi. Considering volume averages 
of the field lead to 



(/i). = -^ 



Ab_ 
E 



3s3CT2(2CTi-fCT2) 



s3(2ai-|-a2)(2CT2 + fTi)-2(5a)2 



(D2) 



where Sa = cj i — CT2 and {fi)x — {fi)y ~ 0. Now expand- 
ing Eqn. (Dl) in powers of Sa gives. 



CTe — 0-2 + i6a)p 



3(72 



[Safp- 



osVi 



{Safp. (D3) 




(a) (b) 

FIG. 19. Periodic cellular models: (a) Sheet-like cell; (b) 
Bond/node-like cell. 

Similarly Brown's formula |5q| to the same order gives, 

a,^a2 + {Sa)p - ^{5afp + ^-^^{5afp. (D4) 
0G2 



9ai 



Equating similar terms leads to CiIp=o = s 



h 



Points representing this result are plotted in Fig. H and 
confirm prior calculations of C,\ . It should also be possi- 
ble to calculate the first order correction, i9Ci/i9p|p=o, by 
calculating (7e to Oij?') [^,|6|. Since Ci is observed to 
have a linear behaviour over a wide range of p pq | (see 
Fig. II) this would provide a good estimate of ^i- Also 
note that r\\ |p=o can be derived using similar methods. 

APPENDIX E: PERIODIC CELL MODELS 

To explicitly demonstrate the effect of pore shape on 
effective conductivity we estimate u^, for several periodic 
networks exhibiting sheet-like, grid-like and node/bond- 
like cells. Consider a structure comprised of periodic rep- 
etitions of the unit cell shown in Fig. ^(a). Defining 
X = hjEl the volume fractions of each phase are given by 
p = 1 - (1 - a:)3 ~ 3x and g = (1 - a;)3 ~ 1 - 3a;. Con- 
sider the behaviour of the model if a\ '^ ai- In this case 
most of the current would flow through the solid faces 
of the cell which are aligned in the direction of current 
flow. The volume fraction of these elements of the cell is 
P\ = 2x — X- . The remaining current would pass through 
a layer of phase 1 (volume fraction p^ = x(\ — x)^) and 
the cell core of phase 2 (volume fraction q). Treating 
each of these mechanisms as conductors in parallel we 
have CTe — piCTi -I- (p2 + 9)0'*, where u* is conductivity of 
the central leakage pathways. Assuming each of the ele- 
ments of these pathways act as conductors in series gives 
CT* = (p2 + q){p2(y\^ + qd^^Y^- This leads to 



{2x 



0-2 



Ao^ 



TPCTl 



{\-xfo2 



x{a2lcri) + {l-x) 



-pa2 IH 

3 V CTi 



where the approximation holds for p <^ 1. Finally, 
(7e ~ Icip in the case cr2 = 0. 
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In a similar way a 'toy' model can be defined to qual- 
itatively demonstrate the effect that necks/throat have 
on the effective conductivity. A cross section of the unit 
cell of a node/bond model is shown in Fig. [l9|(b). The 
central cube has side length x™ and the six arms have a 
square cross section of side length x. Taking the cell to 
have unit width we have: p — x^™ + 3{1 — x™)x'^ (m < 1) 



and q — I — p. If (72 = then most of the current will 
flow through the bonds parallel to the direction of the 
applied field. Therefore, a^ ~ aix^. In the case m — 1 
a uniform grid results and CTe ~ jP'^'i to leading order 
in p. For jn = 1/3 a node/bond geometry results and 
(Te c^p'^ai. 
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